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Abstract 

We  present  transient  and  asymptotic  reliability  indices  for  a  single-unit  system  that  is  subject 
to  Markov-modulated  shocks  and  wear.  The  transient  results  are  derived  from  the  (transform) 
solution  of  an  integro-differential  equation  describing  the  joint  distribution  of  the  cumulative 
degradation  process  and  the  state  of  the  modulating  process.  Additionally,  we  prove  the  asymp¬ 
totic  normality  of  a  properly  centered  and  time-scaled  version  of  the  cumulative  degradation  at 
time  t.  This  asymptotic  result  leads  to  a  simple  normal  approximation  for  a  properly  centered 
and  space-scaled  version  of  the  system’s  lifetime  distribution.  Two  numerical  examples  illustrate 
the  quality  of  the  normal  approximation. 
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1  Introduction 


In  this  paper  we  investigate  transient  and  asymptotic  reliability  indices  for  a  single-unit  system 
that  degrades  over  time  due  to  normal  wear  induced  by  its  operating  environment  and  randomly- 
occurring  shocks  that  cause  additional  damage  to  the  system.  The  wear  rate  and  shock  arrival 
rate  are  both  modulated  by  an  external  process  (the  random  environment)  which  is  modeled  as 
a  continuous-time,  regenerative  stochastic  process  {Zt  :  t  >  0}.  At  time  zero,  the  device  begins 
its  lifetime  in  perfect  working  order  but  degrades  over  time  under  the  influence  of  its  random 
environment  and  shock  damage  until  the  cumulative  degradation  reaches  a  fixed  threshold  x,  at 
which  time  it  fails.  We  denote  this  random  time  until  failure  by  Tx.  Derived  herein  are  transform 
expressions  for  the  cumulative  distribution  function  (c.d.f.)  of  Tx  and  its  nth  moment  (n  >  1)  when 
the  environment  is  assumed  to  evolve  as  a  time- homogeneous  Markov  process.  The  main  transient 
results  are  obtained  by  showing  that  the  joint  distribution  of  the  cumulative  degradation  level 
and  the  environment  state  conditionally  satisfies  an  integro-differential  equation.  This  equation 
leads  to  the  Laplace-Stieltjes  transform  (LST)  of  the  lifetime  distribution  function  and  the  LST 
of  its  moments.  We  additionally  investigate  a  scaled  version  of  the  unit’s  lifetime  with  the  aim 
of  developing  simple  approximations  that  allow  us  to  circumvent  (numerical)  Laplace  transform 
inversion.  The  asymptotic  results  are  obtained  by  analyzing  the  degradation  process  as  a  cumulative 
stochastic  process  which  is  shown  to  obey  a  central  limit  theorem. 

Systems  that  incur  degradation  due  to  their  operating  environment  and  shock  damage  can  be 
identified  in  a  variety  of  settings.  For  instance,  low-observable  technologies  are  commonly  employed 
to  ensure  that  weapons  systems  are  difficult  to  detect,  track  and  engage.  These  “stealth”  systems  are 
usually  coated  with  a  radar- absorbing  material  (RAM)  for  energy  absorption  and  cancelation.  RAM 
coatings  must  maintain  precise  tolerances  for  an  array  of  electrical  properties,  most  importantly  the 
permeability,  and  consequently,  the  near-field  reflectivity.  However,  RAM  is  extremely  sensitive  to 
the  ambient  environment  and  suffers  wear  from  normal  usage  and  exposure  to  the  elements  (e.g., 
sun,  rain,  hail,  etc.).  Additionally,  the  coating  suffers  random  strikes  (e.g.,  scratches  resulting  from 
aerial  refueling,  bird  strikes  and  pebble  strikes  during  take-off  and  landing)  that  cause  damage 
to  the  external  skin,  degrading  the  aircraft’s  near- field  reflectivity.  Assessing  the  impact  of  wear 
and  random  damage  on  the  near- field  reflectivity  of  the  RAM  coating  is  very  difficult,  but  can 
be  done  through  a  detailed  and  involved  inspection  process.  Once  the  degradation  of  the  RAM 
coating  reaches  or  exceeds  a  significant  level,  the  aircraft  must  be  taken  out  of  service  and  repaired. 
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Another  example  of  a  device  that  is  subject  to  a  time-varying  operating  environment  and  random 
damage  is  a  cutting  tool  on  manufacturing  equipment  (e.g.,  the  cutting  tool  on  a  lathe).  The  tool  is 
subject  to  different  wear  rates  that  may  depend  on  several  factors  including  the  composition  of  the 
work  piece,  the  cutting  speed  and  even  variability  among  machine  operators.  The  cutting  tool  may 
also  encounter  defects  in  the  material  that  cause  excessive  tool  wear.  Once  the  tool’s  degradation 
level  reaches  a  critical  threshold,  it  is  unable  to  maintain  engineering  specifications  and  must  be 
replaced. 

Stochastic  failure  models  that  attempt  to  capture  the  impact  of  a  randomly  evolving  envi¬ 
ronment  have  been  examined  extensively  in  the  applied  probability  literature.  An  excellent  survey 
contributed  by  Singpurwalla  [25]  presents  a  number  of  models  with  various  attributes.  Much  earlier 
work,  due  to  Esary,  et  al.  [7]  provided  several  results  for  both  wear  and  shock  processes.  That  work 
has  been  extended  in  numerous  directions  by  several  authors.  For  instance,  Qinlar  [3]  generalized 
most  of  the  models  of  [7]  by  demonstrating  that  the  joint  process  of  the  unit’s  wear  level  and  the 
state  of  its  ambient  environment  may  be  considered  as  a  Markov-additive  process  and  gave  several 
such  examples.  The  first  example  considered  the  case  when  the  random  environment  is  a  general 
Markov  process  with  a  finite  state  space  and  the  wear  is  assumed  to  increase  as  a  Levy  process. 
Additionally,  random  shocks  were  assumed  to  occur  at  environment  transition  epochs.  The  second 
example  (see  [3],  pp.  201-202)  is  similar  to  the  problem  discussed  here  where  the  cumulative  wear 
is  a  continuous,  additive  functional  of  the  operating  environment,  and  the  first  time  to  failure  is  a 
first  passage  time  for  the  degradation  process.  However,  that  model  did  not  include  shocks. 

Though  Rade  [23],  Shanthikumar  and  Surnita  [24],  and  Nakagawa  [21]  provide  extensions  to  the 
models  of  [7]  they  did  not  incorporate  the  effect  of  the  unit’s  operating  environment  on  reliability 
measures.  Ebrahimi  [6]  investigated  survival  functions  for  a  cumulative  damage  shock  model  with 
a  critical  threshold  and  provided  some  stochastic  ordering  results.  Li  and  Luo  [19]  considered  a 
Markov-modulated  shock  process  wherein  the  shock  inter-arrival  times  and  the  random  shock  dam¬ 
age  are  both  governed  by  a  Markov  chain.  They  obtain  reliability  bounds  when  the  inter-arrival 
times  have  heavy-  or  light-tailed  distributions.  Their  degradation  model  does  not  include  a  con¬ 
tinuous  wear  component.  Mallor  and  Orney  [20]  considered  a  generalized  shock  process  and  study 
some  limiting  properties.  Igaki  et  al.  [12],  Skoilakis  [26],  and  Kharoufeh  [13]  present  degradation 
models  that  include  the  influence  of  a  random  environment  but  do  not  consider  random  shocks. 

Klutke  et  al.  [18]  examined  the  availability  of  an  inspected  system  whose  inter-inspection  times 
and  wear  rates  are  random.  Subsequently,  Klutke  and  Yang  [17]  derived  an  availability  result 
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for  a  system  subject  to  constant  linear  degradation,  shocks  and  a  deterministic  inspection  policy. 
Kiessler  et  al.  [16]  investigated  the  limiting  average  availability  of  a  system  whose  time- varying  wear 
rates  are  governed  by  a  continuous-time  Markov  chain.  Kharoufeh  et  al.  [15]  extended  the  model 
of  [16]  by  including  damage-inducing  shocks  that  arrive  according  to  a  time-homogeneous  Poisson 
process  and  deriving  the  Laplace-Stieltjes  transforms  of  a  few  transient  reliability  indices.  However, 
they  did  not  consider  the  asymptotic  behavior  of  these  indices.  Ebrahimi  [5]  considered  a  system 
whose  degradation  is  comprised  of  a  continuous  wear  component  as  well  as  jumps.  The  properties 
of  the  model  were  investigated  and  bounds  were  established  for  the  reliability  function.  However, 
the  random  environment  in  that  model  corresponds  to  a  shock  process  that  is  superimposed  on  a 
gamma  wear  process. 

Our  work  here  extends  the  results  of  [15],  [13]  and  [16]  in  several  ways.  First  we  consider  a 
model  in  which  the  rate  of  continuous  wear  and  the  arrival  rate  of  damage-inducing  shocks  are 
both  modulated  by  an  external  environment.  Initially,  we  focus  on  characterizing  the  transient 
reliability  indices  (i.e.,  when  the  critical  degradation  threshold  is  finite)  in  the  spirit  of  [15,  13]. 
We  provide  a  characterization  of  the  transient  distribution  and  the  lifetime  moments  of  the  single¬ 
unit  system  in  the  form  of  Laplace-Stieltjes  transforms  that  are  amenable  to  numerical  inversion. 
However,  our  aim  is  to  go  beyond  these  transform  solutions  to  gain  further  insight  into  the  behavior 
of  the  degradation  process  (and  the  resulting  lifetime  distribution)  in  an  asymptotic  sense  (as  the 
threshold  goes  to  infinity).  To  this  end,  we  state  and  prove  several  limit  theorems  related  to  a 
time-scaled  version  of  the  degradation  process  and  a  space-scaled  version  of  the  unit’s  random 
lifetime.  Specifically,  we  prove  the  asymptotic  normality  of  these  properly  scaled  quantities.  The 
asymptotic  results  serve  as  simple  approximations  for  the  lifetime  distribution  and  may  be  useful  for 
degradation-based  reliability  models  such  as  those  described  by  Gebraeel  et  al.  [9]  and  Kharoufeh 
and  Cox  [14]  among  others. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  provides  the  model  description 
and  the  notation  used  throughout  the  paper.  In  section  3,  we  characterize  the  lifetime  distribution 
and  the  nth  moment  (n  >  1)  of  the  lifetime  in  a  transient  sense  (i.e.,  when  the  degradation 
threshold  x  is  finite).  These  characterizations  are  in  the  form  of  Laplace-Stieltjes  transforms. 
Section  4  investigates  a  time-scaled  version  of  the  degradation  process  and  a  space-scaled  version 
of  the  system’s  lifetime,  both  of  which  are  shown  to  obey  a  central  limit  theorem.  Two  illustrative 
examples  are  included  in  section  5,  while  section  6  provides  a  few  concluding  remarks. 
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2  Model  Description 


In  this  section,  we  describe  our  mathematical  model  for  a  single-unit  system  subject  to  environment- 
induced  wear  and  shocks.  All  random  variables  are  defined  on  a  common  and  complete  probability 
space  (D,A,P).  A  single- unit  system  is  placed  into  service  at  time  zero  in  perfect  working  order. 
The  system  accumulates  degradation  until  a  deterministic,  critical  threshold  value  x  is  reached  or 
exceeded,  at  which  time  the  system  is  said  to  be  failed.  Denote  this  random  first  passage  time  by 
Tx.  We  assume  throughout  that  Tx  is  proper,  i.e. ,  for  each  x  >  0, 


lim  P(TZ  <  t)  =  1. 

t— XX) 

In  fact,  the  dynamics  of  our  model  ensure  that  Tx  is  bounded.  The  degradation  accrued  over  time 
by  the  system  is  attributed  to  environment-induced  wear,  as  well  as  shocks  that  occur  at  random 
time  intervals.  Both  the  wear  and  shock  arrival  rates  are  modulated  by  an  external  stochastic 
process  (the  random  environment).  The  wear-  and  shock-inducing  mechanisms  are  assumed  to 
behave  independently  of  one  another;  however,  they  share  a  common  dependence  on  the  state  of 
the  random  environment. 

The  random  environment  is  assumed  to  be  a  regenerative  process,  namely  an  irreducible 
continuous-time  Markov  chain  (CTMC),  Z  =  {Zt  :  f  >  0},  on  a  finite  state  space  S  =  {1,  2, . . .  ,£}. 
The  CTMC  has  an  infinitesimal  generator  matrix  Q,  initial  distribution  vector  a  and  transition 
probability  functions  Trt.j(t),  i,j  £  S  which  comprise  the  matrix  11(f)  =  [7r,j(t)].  The  limiting  dis¬ 
tribution  of  Z  is  7T.  Define  positive  functions  r,  A  :  S  — >  R+  such  that  whenever  Zt  =  i,  the  system 
wears  linearly  at  a  rate  r(i),  r(i)  >  0,  and  shocks  occur  according  to  a  Poisson  process  at  rate  A*, 
Aj  >  0  for  i  =  1,2, . . .  That  is,  both  the  wear  rates  and  shock  arrival  rates  are  Markov- modulated. 
For  notational  convenience,  let  =  diag(r(l),  r(2), . . . ,  r(£ ))  and  A  =  diag(Ai,  A2, . . . ,  A^).  While 
the  assumption  of  linear  wear  appears,  on  first  glance,  to  be  restrictive,  non-linear  wear  paths  can 
be  effectively  approximated  by  piece- wise  linear  paths  as  demonstrated  by  Kharoufeh  and  Cox  [14] . 

Next  we  describe  the  total  degradation  process.  The  cumulative  wear  up  to  time  f,  denoted  by 


Wt,  is 


Wt  =  W0  + 


r(Zu)du 


(1) 


where  we  assume  Wo  =  0,  and 

r(Zu)\du  <  00  a.s. 

so  that  Wt  is  well  defined  for  each  f  >  0.  The  process  W  =  {Wt  :  f  >  0},  is  a  cumulative  stochastic 
process  (i.e.,  an  additive  functional  of  a  regenerative  process)  as  defined  in  [2,  11,  10,  27].  The 
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system  is  also  damaged  by  shocks  that  arrive  according  to  a  Poisson  process  with  rate  A  ;  whenever 
Zj  =  i.  The  damage  caused  by  an  individual  shock  is  assumed  to  be  relatively  small;  however, 
the  cumulative  effect  of  small  shocks  may  be  significant,  as  in  the  case  of  fatigue  deterioration 
resulting  from  mechanical  vibrations.  Denote  by  Nf  the  number  of  shocks  occurring  up  to  time  t. 
The  corresponding  counting  process  {Nt  :  t  >  0}  is  a  Markov-modulated  Poisson  process  (MMPP) 
with  £  Poisson  arrival  rates,  Ai,  A2, . . . ,  A  £  (see  [8]  for  further  details).  The  damage  caused  by  the 
nth  shock  is  a  (nonnegative)  random  variable  Yn,  and  {Yn{^=l  is  an  i.i.d.  sequence  with  proper 
c.d.f.  Fy(y)  =  P(Y"  <  y).  mean  fj,  =  E(Yi)  and  variance  o\-  =  Var(Yi).  We  assume  throughout  that 
E(Yi)  <  00  and  E(Y12)  <  00.  The  cumulative  damage  due  to  shocks  up  to  time  t  is  a  nonnegative 
random  variable, 

Nt 

Pt  =  J2Yn’  t^°-  (2) 

n=0 

The  total  degradation  accrued  by  the  system  up  to  time  t  is  the  sum  of  degradation  due  to  wear 
and  that  due  to  shocks  given  by 

rt  Nt 

Xt=  r(Zu)du  +  y^jYn,  t>  0.  (3) 

Jo  n= 0 

The  positivity  of  the  degradation  rates,  r(l),r(2), . . .  ,r(£),  ensures  that  the  sample  paths  of  X  = 
{Xt  :  t  >  0}  are  monotonically  increasing  almost  surely  (a.s.),  and  consequently,  that  events 
{Xt  <  x}  and  {Tx  >  t}  are  equivalent.  The  system’s  random  lifetime  is  given  by 

Tx  =  inf {t  >  0  :  Xt  >  x},  (4) 

or  the  first  time  the  degradation  process  X  reaches  or  exceeds  x.  Let  G(x,  t)  =  P (Tx  <  t)  = 
1  —  P(W  <  x )  denote  the  unconditional  c.d.f.  of  the  unit’s  lifetime,  and  let  its  nth  moment  be 
denoted  by  E(T”)  for  n  >  1.  This  paper  is  concerned  with  two  primary  aims.  First,  we  characterize 
the  transient  versions  of  G(x,  t)  and  E(T"),  n  >  1  and  their  conditional  counterparts.  By  transient, 
we  mean  for  all  x  such  that  0  <  x  <  00.  Second,  we  investigate  a  centered  and  time-scaled  version 
of  Xt  as  t  — ►  00,  and  a  centered  and  space-scaled  version  of  Tx  as  x  — ►  00,  to  construct  asymptotic 
approximations  for  the  transient  indices.  Sections  3  and  4  provide  our  main  results,  and  section  5 
presents  a  few  numerical  illustrations. 

3  Transient  Analysis 

This  section  provides  the  transient  reliability  indices  for  a  single- unit  system  whose  degradation 
evolution  is  described  by  a  Markov-modulated  shock  and  wear  process.  Provided  are  expressions 
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for  the  system’s  lifetime  c.d.f.,  as  well  as  the  moments  of  the  lifetime,  in  the  form  of  Laplace- 
Stieltjes  transforms  (LSTs).  Let  us  first  introduce  the  following  definitions  and  notation.  The 
complementary  c.d.f.  of  the  lifetime  is 

R(x,t)  =  P(Xf  <  x)  =  1  —  G(x,t).  (5) 

Define  the  joint  probability  distributions 

Rij(x,t)  =  F(Xt  <  x,Zt  =  j\Z0  =  i),  i,jeS.  (6) 

Our  first  main  result  characterizes  the  joint  distribution  of  the  process  (T,  Z),  conditioned  on  the 
initial  state  of  the  environment.  This  result  intuitively  extends  Theorem  1  of  [15]  to  the  case  of  a 
time-varying  (modulated)  shock  arrival  rate. 


Theorem  1  For  each  i,j  €  S,  the  distribution  function  Rjj(x,t)  verifies  the  partial  integro- 
differential  equation 


dRij(x,t )  dRij(x,t ) 


dt 


+ 


dx 


r(j )  =  A j  ([ Rij(-,t )  *  TV]  (x)  -  Ri,j(x,  t))  +  Y  Qk,jRi,k(x’ f)  ( 7 ) 


k= l 


for  x  >  0,  t  >  0  where  (*)  denotes  the  convolution  operator. 


Proof.  The  proof  is  similar  to  that  of  Theorem  1  in  [15].  Considering  only  the  cumulative 
wear  of  the  system  first,  let  Vij(x,t)  =  P (W*  <  x,  Zt  =  j\Zo  =  i )  and  let  e  >  0  denote  a  small 
time  increment.  Because  Z  is  a  temporally  homogeneous  Markov  process  and  is  independent  of 
the  degradation  process  T,  we  may  write 


Vij(x,t  +  e) 


P (Wt+e  <  x ,  Zt+€  =  j\Z0  =  i ) 

l 

Y.  p (Wt+e  <  X,  Zt+e  =  j\Zt  =  k,  Z0  =  i)  P (Zt  =  k) 

k= l 

l 

Y^izt+e  =  j\Zt  =  k,Z0  =  i)  P (Wt+e  <x\ Zt  =  k,  Z0  =  i )  P (Zt  =  k ) 

k=  1 
l 

]TP(Zt+e  =  j\Zt  =  k )  P(Wt+e  <x,Zt  =  k\Z0  =  i) 

k= 1 

l 

Y,  VTjtj(e)  P {Wt+e  <x,zt  =  k\Z0  =  i) 

k=  1 
l 

Y71 'k,j{t)Vi:k(x  -  r{k)e,t).  (8) 

k= 1 
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Next  we  incorporate  the  impact  of  Markov-modulated  Poisson  shocks  occurring  at  rate  A zt-  Let 
Nf  denote  the  random  number  of  shocks  up  to  time  t.  The  transition  probability  functions  for  the 
Z  process,  7 Tij(e),  i,j  €  S,  can  be  written  as 

7 =  Sij  +  eqitj  +  o(e)  (9) 


where  5ij  assumes  the  value  1  when  i  =  j  and  0  when  i  ^  j.  It  is  well  known  that 


P (Nt+e  -Nt  =  0\Zt  =  k) 

—  1  —  Afce  +  o(e), 

(10a) 

P  (Nt+e  —  Nt  =  l\Zt  =  k) 

=  A  fce  +  o(e), 

(10b) 

F(Nt+e  —  Nt>  2\Zt  =  k) 

=  o(e). 

(10c) 

The  magnitude  of  cumulative  damage  caused  by  n-independent  shocks  is  given  by 

n 

Pn  =  J2Yi- 

i= 1 

Since  the  shock  magnitudes  form  an  i.i.d.  sequence  with  distribution  function  Fy,  we  note  that 

nPn  <y)  =  F<?\y), 

(n) 

where  Fy  denotes  the  n-fold  convolution  of  Fy  with  itself.  Conditioning  on  the  number  of  shocks 
in  the  interval  ( t ,  t  +  e)  and  the  magnitude  of  damage  due  to  those  shocks,  we  may  now  write 

i  OO 

Rij(x,t  +  e)  =  Y:  ni,k(e)  F( Nt+e  -  Nt  =  n\zt  =  k ) 

k= 1  n= 0 

x  f°°  ¥(Xt+e  <  x,  Zi+e  =  j\Zt  =  k,Y  =  y)Fp\dy). 

Jo 


(11) 


/  0 

Substituting  equations  (8)  and  (10)  into  (11)  we  have 

l 


Ri,j{x,  t  +  e)  =  Y  x 

k= 1 

^ Ri,k(x  -  r(k)e,t)(  1  -  Xke)  +  A ke  J  Ri,k{x  -  r(k)e  -  y,t)FY(dy)Sj  +  o(e).  (12) 
Using  (9)  in  (12)  and  simplifying  gives 

l 

Ri,j{x,  t  +  e)  =  (1  -  A ke)Ritj(x  -  r(j)e ,  t)  +  e  ^(1  -  A ke)qk,jRi,k{x  ~  r(k)e,  t ) 

k= 1 

poo 

+  Aie  /  -  r(j)e  -  y,i)iV(dy) 

Jo 

I  poo 

+  e2  Y  XkQk,j  /  -  r(fc)e  -  y,t)FY(dy)  +  o(e).  (13) 

fe=i  Jo 


Rearranging  and  simplifying  the  terms  of  (13),  dividing  by  the  time  increment  e,  and  letting  e  j  0, 
shows  that 


dRi,j(x,t)  dRij(x,t) 

—a^  +  —a^rb)  = 

I  /*°° 

-  A jRij(x,  t)  +  y~]  Qk,jRi,k(x,  t)  +  A j  /  Ri,j{x  -  y,  t)FY(dy)  (14) 

k= i  Jo 

where  the  right-most  term  of  (14)  is  the  convolution  of  the  distributions  Rij  and  Fy. 


The  differential  equation  (7)  describes  the  joint  (spatial)  evolution  of  the  degradation  process  with 
the  (temporal)  evolution  of  the  unit’s  operating  environment.  The  impact  of  shocks  is  captured  by 
the  convolution  term  on  the  right-hand  side  of  (7). 

The  system  of  equations  can  be  written  in  matrix  form  to  make  the  solution  procedure  more 
transparent.  Define  the  i  x  i  matrix  R(x,  f)  =  [Rjj(x,  f)]  and  recall  that  A  =  diag(Ai,  A2, . . . ,  A^). 
In  matrix  notation,  (7)  may  be  written  as 

+  9KQXx,t]Rd  =  ([R M)  *  Fy]  (x)  -  R(x,  t))  A  +  R(x,  t) Q  (15) 

where  [R(-,t)  *  TV]  (x)  denotes  the  convolution  of  each  entry  of  R(x,f)  with  Fy.  The  initial  prob¬ 
ability  vector  of  the  environment  process  Z  is  ol  =  [cq]  where  oii  =  P(^o  =  i),  i  €.  S.  Let 
Pj(Al)  =  W{A\Zq  =  i )  for  A  e  A,  and  let  Ej(T")  =  E(T”|Zo  =  i).  The  vector  e  will  denote  a  column 
vector  of  ones,  and  e,;  is  a  column  vector  whose  ith  entry  is  unity  and  all  others  are  zero.  The 
unconditional  c.d.f.  of  the  system’s  lifetime  is  given  by 

G(x,t )  =  P (Tx  <  t)  =  1  —  aR(x,t)e,  (16) 


and  the  conditional  c.d.f.  is  given  by 


Gi(x,t )  =  P i(Tx  <t)  =  1  -  e'R(x,  t)e 


(17) 


where  e\  denotes  the  transpose  of  e*.  Let  the  Laplace- Stieltjes  transforms  of  G(x,t)  and  Gi(x,t) 
with  respect  to  x  be 


and 


G(u,t)  =  /  e  uxGY(dx,  t).  Re(u)  >  0, 
Jo 


Gi(u,t)  =  /  e  wxGi(dx,  t),  Re(it)  >  0, 
Jo 
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respectively.  Define  an  l  x  i  diagonal  matrix,  Fd(u),  with  each  diagonal  entry  identically  equal  to 

FY(u)  =  /  e~uyFY(dy),  (18) 

Jo 

the  LST  of  Fy  with  respect  to  y.  The  following  result  generalizes  the  lifetime  distributions  of  [16], 
[13]  and  [15]  to  the  case  of  a  Markov-modulated  shock  and  wear  process. 

Theorem  2  The  Laplace- Stieltjes  transforms  of  the  c.d.f.  ’s  G(x,t)  and  Gi(x,t),  with  respect  to  x 
are,  respectively, 

G(u,t)  =  1  —  cx  exp  (X('u)t)  e  (19) 

and 

Gi(u,t )  =  1  —  e(  exp  (X(it)f)  e,  (20) 

where  X(u)  =  Q  +  (F d(u)  -  I)A  -  uRd,  I  is  the  identity  matrix  and  Re(n)  >  0. 

Proof.  Define  the  matrix  LST  of  R(a;,  t)  with  respect  to  x  as 

~  r 00 

R  (u,t)=  /  e~uxR(dx,t).  (21) 

Jo 

Taking  the  LST  of  both  sides  of  (15)  with  respect  to  x  yields  a  first  order,  ordinary  differential 
equation  in  t, 

dR^  +  R (u,  t)  [uRd  +  (i  -  Fd(«))  A  -  Q]  =  0.  (22) 

Applying  the  initial  condition,  R(it,  0)  =  I  and  rearranging  terms,  this  ordinary  differential  equation 
has  the  obvious  solution, 

R(u,  f)  =  exp  (^Q  —  (I  —  F d{u)) X  —  uRd^j  t  . 

Thus,  the  LST  of  G(x,  t )  is  given  by 

G(u,  t ) 

The  LST  of  the  conditional  c.d.f.,  Gi(x,t),  is  obtained  by  replacing  a  by  e\.  ■ 

We  next  examine  the  unconditional  and  conditional  moments  of  the  system  lifetime.  To  this 
end,  let  Tf  be  the  nth  power  of  the  random  variable  Tx,  and  respectively  denote  the  unconditional 
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and  conditional  moments  of  Tx  by  E (T”)  and  E i(T£),  n  >  1.  Define  the  LST  of  E(T”)  with  respect 
to  x  by 

E(T”)  =  /  e— dE(T”), 

Jo 

and  the  LST  of  E,;(T”)  as 

^  r°° 

E  i(2£)=  /  e_ua;dE(Tx|Zo  =  i). 

Jo 

Corollary  1  For  eac/i  n  >  1,  t/ie  Laplace- Stieltjes  transform  o/E(T”)  with  respect  to  x  is  given  by 

E{TZ)  =  n\-a(uRd  +  {l-FY{u))\-Qyne,  (23) 

and 

~  ~  \  —  n 

Ei(TZ)  =  n\-e'i(uRd+(l-FY(u))\-Q)  e  (24) 

with  Re(rt)  >  0. 


Proof.  The  result  can  be  shown  by  first  defining  the  Laplace-Stieltjes  transform  of  G(u,  t ) 
(with  respect  to  t)  given  by 

roc  ~ 

G(u,s)=  /  e~stG(u,dt). 

Jo 

Apply  standard  matrix  calculus  operations  (cf.  Neuts  [22]),  we  obtain  the  first  result  by 


Era  =  (-1) 


tdG(u,s) 


d  sr 


5—0 


The  conditional  version  (24)  is  obtained  by  simply  replacing  G(u,s)  by  Gi(u,s )  in  the  above 
expression.  ■ 

We  note  that  if,  for  some  constant  A  (A  >  0),  A,  =  A  for  all  i  G  S,  the  results  of  Theorem  2 
and  Corollary  1  reduce  to  the  corresponding  LSTs  of  [15].  Furthermore,  if  A i  =  0  for  all  i  G  S, 
the  model  reduces  to  a  Markov-modulated  wear  process  with  no  shocks  as  described  in  [13]  and 
[16].  The  transient  results  presented  here  are  amenable  to  numerical  Laplace  transform  inversion; 
however,  numerical  inversion  can  be  unstable  when  x  is  either  very  small  or  very  large,  or  when  \S\ 
is  large.  Alternatively,  we  may  consider  the  limiting  behavior  of  the  degradation  process  with  the 
aim  of  constructing  asymptotic  distribution  and  moment  approximations. 
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4  Limiting  Behavior 


We  now  examine  the  cumulative  degradation  process  and  lifetime  distribution  in  the  asymptotic 
regime  (i.e.,  when  t  — >  oo  and  x  — >  oo,  respectively).  As  before,  the  row  vector,  7r  =  [7^]  for 
i  =  1,2, ...  ,£,  denotes  the  stationary  distribution  of  the  environment  process  Z.  Let 

A  =  (R-d  +  /xA)e, 

and  note  that  7tA  is  the  long-run  average  rate  of  degradation  (including  wear  and  damage  from 
shocks).  Intuitively,  we  would  expect  that,  when  x  is  large,  the  average  time  for  the  degradation 
process  to  first  reach  level  x  should  be  very  close  to  x/tvA.  Proposition  1  confirms  this  intuition 
by  establishing  convergence  in  the  mean  of  the  scaled  r.v.  x~lTx  as  x  — >  00. 


Proposition  1  As  x  — >  00, 


E(TX 


1 


X  7tA 

Proof.  To  prove  this  result,  we  first  note  that  as  u  — >  0 


(25) 


Fy(u )  —  1  d  ~ 

- >  Fy  «) 

u  du 


=  ~fi, 


u= 0 


so  that 


1  —  Fy(u) 


u 


=  n  +  e(u) 


where  ue(u)  =  o(u).  Therefore,  using  (23), 


lim  a[uRd  +  (1  —  Fy{u)) A  —  Q] 

u — >0 


1  =  lim  u 

u — ^0 


f  1-Fy(u)\\  n 

u  Kd  +  - — -  A  -  Q 


1  -1 


u 


=  lim  u  [u(Rrf  +  //A)  —  Q  +  o  (u)A] 


-1 


it— >0 


^lim  u  1  [u(Rrf  +  /xA)  —  Q  +  o(xx)A]  j 


-1 


=  (  lim  u  1  [u(Rd  +  /xA)  —  Q 

\u^  0 


=  lim  u  [u(Rrf  +  /xA)  —  Q] 

u—>  0 


-1 


-1 


Next  we  make  use  of  the  asymptotic  properties  of  the  Laplace  transform  to  obtain, 

®(tx 


lim 


IE— >00  X 


=  lim  uE(Tu) 

u— >0 


=  lim  «a[«Rd  +  (1  —  Fy(u))X  —  Q]  1e 

u—>  0 

=  lim  aji[«(Rd  +  fiX)  —  Q]_1e 

u— >  0 


-1, 

— — 1/ 


=  lim  olu[u1  —  (Rd  +  yU A)  Q]  (R<j  +  jiX)  ie. 

u— >  0 


12 


Since  (R^  +  pX)  1  is  diagonal,  Q  =  (R^  +  pX)  XQ  is  a  generator  matrix  for  a  CTMC  {Zt  :  t  >  0} 
with  transition  probability  functions  11(f)  and  limiting  distribution, 

tt(R  d  +  //A) 


7T  = 


7tA 


(26) 


since  (R^  +  pX)  is  positive  recurrent.  The  transition  probability  functions  have  LST 


'L(n)  =  /  e~utdn(t)  =  u(ul  -  Q) 
Jo 

Applying  the  asymptotic  properties  of  the  LST  gives, 


-l 


(27) 


lim  'L(w)  =  lim  fl(f)  =  fl(oo) 

t— >oo 


Therefore,  we  may  write 


lim 


E(TS 


x— >oo  X 


=  lim  cc^(it)(Rrf  +  //A)  1e 

u—>  0 

=  ctn(oo)(Rrf  +  pX)~1e 
1 


7tA 


Before  generalizing  Proposition  1  to  the  nth  moment  (n  >  2),  we  need  the  following  lemma. 

Lemma  1  Let  /(•)  be  a  function  of  exponential  order  on  the  positive  real  line  such  that  f^m\t) 
oo  as  t  — >  oo  for  m  =  0, 1, . . .  ,n  —  1,  and  let  f(s)  denote  its  Laplace- Stieltjes  transform.  Then 

m 


Proof.  It  is  well  known  that 

(d  n+1f\ 


lim  sn  f(s)  =  n!  lim 

v  '  t— oo  tn 


71+1 


c 


-s+r  (»)  =  s”+7»  -  V  s*-1/<»-*+1)(o) 

^  '  fc=l 


where  f*(s)  denotes  the  Laplace  transform  of  /.  Letting  s  — >  0  on  both  sides  of  the  above  expression, 
we  have 

lim  »"+7*(a)  -  /<">«»  =  fcC(^)(*) 

POO 

=  lim  /  e~stf{n+1\t)dt 
s^°  J 0 

=  lim  f1  /(n+1)(t)df 
0.^00  J0 

=  lim  /H(o)-/(")(0). 
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That  is 


lim  snf(.s)  =  limsn+1/*(s)  =  lim  f^n\t)  =  n\  lim 

s^>0  s—>0  t—>oo  t— >oo  tn 

where  the  last  equality  is  obtained  by  n  applications  of  l’Hospital’s  rule. 


Proposition  2  For  each  n  >  2,  as  x  — >  oo, 

E(T^ 


(7rA)n 


Proof.  Lemma  1  gives  us 


lim  unE(T™)  =  n\  lim  ® 
u— >0  x— >oo  xn 


(28) 


That  is, 


E(Tn)  1 

lim  =  -r  lim  unE(T”) 


x— >oo  x 


n! 

=  lim  a  Pu(ul  —  (R^  +  ^A)_1Q(Rrf  +  /iA)_1)l n  e 

L  J 


=  a 


n(oo)(Rrf  +  /xA)  1 
1 


(7rA)r 


Proposition  2  provides  some  insight  as  to  the  behavior  of  x  1TX  as  the  degradation  threshold  x 
gets  large.  In  particular,  we  see  that  as  x  — >  oo, 


However,  a  stronger  result  can  be  proved.  Note  that  the  environment  process  {Zt  :  t  >  0}  is 
regenerative  with  regeneration  epochs  Sq  =  0  and 


Sk+ 1  =  nrin{t  >  Sk:  Zt  =  Z0l  Zt-  /  Z0},  k  >  0. 

Let  f  =  7rRrfe  denote  the  long-run  average  wear  rate,  let  k  =  l/E(S'i)  be  the  average  cycle  frequency 
and  A  =  7rAe  is  the  long-run  average  shock  arrival  rate.  Using  these  quantities,  we  now  prove  that 
x~lTx  converges  almost  surely  to  a  fixed  value. 


Theorem  3  As  x  — ■»  oo, 


Tx  1 

X  7T  A 


a.s. 


(29) 
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Proof.  Let  Uq  =  Ws0  =  0  and  Uk+\  =  Wsk+l  —  1 Vsk  f°r  k  >  0.  We  note  that  {Uk}^=i  is 
an  i.i.d.  sequence  of  random  variables.  Thus,  by  the  strong  law  of  large  numbers  (SLLN)  and  the 
renewal  reward  theorem, 


^  =  W  •  i&(WSl)  =  r  a.s. 

On  On  Tl 

k=  1 

We  see  also  by  the  SLLN  that 


Pt 

t 


Nt 

t 


1 

Nt 


E5‘ 


k=l 


A  n 


a.s. 


Combining  the  shock  and  wear  components,  we  have 


r  Xt  ..  Wt  +  Pt 

hm  —  =  hm  - 

t—>  oo  t  t — »oo  t 


lim 

n^oo 


+  lim 

t—>  OO 


Pt 

t, 


r  +  A/U  a.s. 

7tA  a.s. 


Now  note  that  Xtx  =  x  <  Xtx+ i,  and  so 


lim  inf  —  =  lim  =  7tA  =  lim  — >  lim  sup  —  a.s. 

ai  kx)  Tx  x  >0°  Tx  x  >OC)  ±x  x — >oo  Tx 


Therefore, 

X 

lim  —  =  7r  A  a.s., 

x^oo 

which  implies  our  result.  ■ 

Next  our  aim  is  to  characterize  centered  and  scaled  versions  of  Xt  and  Tx  as  t  — >  oo  and  x  — ►  oo, 
respectively.  To  this  end,  define  the  random  variable 

fSl 

ipi(r)=  /  [r(zu)-r]du  =  WSl-rS1, 

Jo 

and  let  a'^v  =  Var(-0i (f))  denote  the  cyclic  wear  variance.  In  the  results  that  follow,  JV(0, 1)  denotes 
a  standard  normal  random  variable  (i.e. ,  one  with  mean  0  and  variance  1 ) ,  and  =>  means  convergence 
in  distribution.  Before  stating  our  main  results,  it  will  be  helpful  to  review  an  important  result  due 
to  Glynn  and  Whitt  [11]  which  applies  to  both  the  shock  and  wear  processes. 


Theorem  4  (Glynn  and  Whitt  [11]).  There  exist  constants  a  and  b  such  that 

t~1/2(Wt  -  at)  =>  bN( 0, 1)  as  t—>  oo 
if  and  only  if ’E(ifi{a))  =  0  and  Var('0i(a))  =  E(£i)62. 
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This  result  establishes  necessary  and  sufficient  moment  conditions  for  the  shock  and  wear  processes 
to  converge  to  a  standard  normal  variable  when  properly  scaled. 


Proposition  3  As  t  — ►  oo, 


where  a 


t-1/2{Xt  -  tv  At)  =>  aN{  0, 1), 


Proof.  To  prove  this  result,  we  can  consider  the  damage  due  to  shocks  and  the  wear  degrada¬ 
tion  separately.  For  the  shock  damage,  note  that  t~lNt  — >  A  a.s.,  so  we  may  conclude  by  Theorem 
2.1  of  [11]  that 

Nt 

t~1/2(Pt  -  liXt)  « r1'2  J2(Yn  -»)=>  ^(0,  A4)  (30) 

n= 0 

as  t  — >  oo .  The  environment  process  Z  is  regenerative  with  the  first  regeneration  epoch  having 
duration  S\.  Therefore,  the  renewal  reward  theorem  shows  that 


E(V’i(r))  =  E 


=  0. 


Furthermore,  we  have  that 

Var(^i(f))  =  k4E(Si). 


Because  V’(r)  satisfies  the  moment  conditions  outlined  in  [11],  we  conclude  that 


t  1/2(Wt  -  rt)  =>  N(0,  na2v) 


(31) 


as  t  oo.  The  result  is  finally  obtained  by  combining  equations  (30)  and  (31)  and  noting  that 
7tA  =  f  +  /iA.  ■ 


Proposition  3  provides  the  means  by  which  to  analyze  a  centered  and  space-scaled  version  of  the 
lifetime  distribution  in  the  asymptotic  regime. 


Theorem  5  As  x  — ►  oo, 


where  a 


a 

7tA 


N (0,1). 
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Proof.  The  equivalence  of  events  {Tx  >  t}  and  {Xt  <  a:}  shows  that 

.  /  X;  —  (nA.)t  x  —  (tvA  )t\ 

P (Tx  >t)=  P(Xt  <  x)  =  P  — - . 

V  G\/t  G\Jt  J 

Let  t,x  — >  oo  in  a  nontrivial  way,  i.e.,  such  that 

x  —  (7tA  )t 


<7 


where  c  is  a  real  constant.  By  Theorem  3,  this  equivalently  means  that 

t  —  xj  (7tA) 


i J a2x/ (7tA)3 


— c, 


as  t,  x  — >  oo.  Since 


P(TX  >  t)  = 


-  x/ (tv A)  t-x/ {n A) 


\J  cy^x  /  (7tA)3  A/cr2a:/(7rA)3  y 


we  have 


lim  P 


Tx  -  xj  (tv  A)  t-x/ (7tA) 


=  lim  P 


X*  —  (7tA  )t  X  —  (7tA  )t 


t.,x^o o  ^  y/ cr2x/ (7tA)3  cr2x/(7rA)3  y  t,x^oo  y  cr  \/i 

Denote  by  $(•)  the  standard  normal  distribution  function.  Then  we  may  write, 


<7 


Vi 


lim  P 


Tx  -  x / (tv A) 


x  *oo  y  y/ cr2x/ (7rA)3 

by  Proposition  3.  Thus, 

Tx  -  x/(tvA 


>  — c  =  lim  P 

/  t—>  OO 


Xt  -  (7tA )t 


a 


Vi 


<  c  =  T(c) 


lim  P 


x^oo  \  y/ <72x/(7tA)3 


<  c  =  1  —  <h(— c)  =  $(c), 


which  implies  our  result. 


Theorem  5  provides  a  simple  normal  approximation  to  the  lifetime  distribution  that  may  be 
used  as  a  surrogate  for  equations  (19)  and  (20).  This  approximation,  along  with  the  asymptotic 
moment  results,  can  be  used  in  a  variety  of  ways.  For  instance,  the  standard  reliability  indices 
(e.g.,  mean  time-to-failure  or  the  reliability  function)  can  be  compared  for  identical,  single-unit 
systems  that  operate  in  different  environments.  Moreover,  the  asymptotic  normality  of  the  lifetime 
distribution  can  be  used  in  the  context  of  sequential  decision  making,  particularly  for  adaptive 
maintenance  planning.  By  inspecting  and  updating  the  unit’s  degradation  level  periodically,  one 
can  prescribe  dynamic  optimal  maintenance  actions  that  evolve  temporally  with  the  unit’s  remain¬ 
ing  useful  lifetime  distribution.  The  remaining  lifetime  distribution  can  be  easily  computed  (via 
look-up  tables)  using  the  asymptotic  normal  distribution,  whereas  it  is  very  cumbersome  to  repeat¬ 
edly  compute  its  transient  counterpart.  The  next  section  illustrates  the  quality  of  the  asymptotic 
approximations  through  two  numerical  examples. 
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5  Numerical  Illustrations 


This  section  provides  two  numerical  examples  that  illustrate  the  usefulness  of  the  simple  asymp¬ 
totic  approximations  of  section  4.  To  this  end,  let  7  =  1/7tA,  and  define  the  random  variable 

Tr  =  (7x)-1/2  (Tx  -  'yx) . 

Theorem  4  asserts  that  Tx  =>•  'yaN(0, 1)  as  x  — >  00.  Define  the  (conditional)  c.d.f.  of  Tx,  given 
that  the  initial  state  of  the  environment  is  i  G  S,  as 

Fi(x,t;  7)  =  P(TX  <  t\Z0  =  i) 

=  P((7x)-1/2(Tx-7x)  <t\Z0  =  i) 

=  Gi(x,ty/ 7T  +  7X).  (32) 

Next  denote  by  <L(i;7,  a)  the  c.d.f.  of  a  normal  random  variable  with  mean  zero  and  variance 
(ycr)2.  Furthermore,  let  T  C  M+  =  [0,oo),  and  denote  the  maximum  absolute  deviation  between 
Fi(x,t\i)  and  <f>(t;7,  a)  by 

sup|Fj(x,t;7)  -  6(t;7,cr)|,  x  >  0. 
ter 

By  Theorem  5,  we  can  assert  that  sup|.Fj(x,  f;  7)  — <h(t;  7,  cr)  |  — >  0  as  x  — >  00.  However,  the  analytical 

ter 

rate  of  convergence  is  unavailable  because  the  distribution  function  Fi(x,t;  7)  can  only  be  obtained 
numerically  via  the  inverse  Laplace  transform, 

Fi(x,  t]  7)  =  CT1  [u~xGi{u,t^pyx  +  7a;))  (33) 

where  T”1  denotes  the  inverse  Laplace  transform  operator.  This  inverse  transform  will  serve  as  the 
baseline  distribution  against  which  we  will  compare  the  asymptotic  normal  approximation  4>(t;  7,  a) 
on  the  fixed  set  T  for  several  (increasing)  degradation  thresholds. 

Additionally,  we  would  like  to  assess  the  quality  of  the  asymptotic  moment  approximations  by 
comparing  them  to  the  transient  moments  given  by 

E (7T)  =  C'1  (u-1  ECO)  ,  n  >  1.  (34) 

The  asymptotic  approximation  of  the  nth  moment  of  Tx  is  given  by 

n 

E(T")f*(^’  (35) 

The  transient  c.d.f.  values  of  (33),  and  the  transient  moments  of  (34),  were  obtained  by  coding 
the  inversion  algorithms  of  Abate  and  Whitt  [1]  in  the  MATLAB  computing  environment. 
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5.1  An  Alternating  Environment 

The  first  example  is  a  simple  environment  which  alternates  between  two  distinct  states  ( S  = 
{1,2}).  Thus,  the  environment  may  be  viewed  as  an  alternating  renewal  process  (or  up-down 
machine)  wherein  the  up  and  down  times  are  mutually  independent.  The  (positive  recurrent) 
generator  matrix  of  the  process  is 

-8.0  8.0 

Q  = 

6.0  -6.0 

and  its  stationary  distribution  is 

tt=[  0.4286  0.5714]- 

The  diagonal  matrix  of  wear  rates  is  =  diag  (1.0, 0.25),  the  diagonal  matrix  of  Poisson  shock 
arrival  rates  is  A  =  diag  (0.25, 2.0),  and  the  vector  A  is  given  by 

1.0625 

A  = 

0.7500 

The  Laplace-Stieltjes  transform  of  the  shock-magnitude  distribution  function  is 

*W«)  =  Re(“)  >  °» 

i.e. ,  the  damage  induced  by  each  shock  is  an  exponential  random  variable  with  rate  parameter  4.0. 
Using  these  input  parameters,  and  assuming  the  environment  is  initially  in  state  1  with  probability 
1.0,  we  can  easily  compute  the  following: 

7=  w  1.13131;  A  =  vrAe  «  1.25000; 

7T  A 

crjV  ~  0.02778;  =  0.06250; 

k  =  l/E(Si)  ~  3.4286;  a  =  \J  na^v  +  A  cry  «  0.41637. 

First  we  compare  the  lower  (first  and  second)  moments  of  Tx  using  equations  (34)  and  (35).  The 
results  are  summarized  in  Table  1. 
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Table  1:  E(T™  |Zq  =  1),  n  =  1,2,  when  Z  has  two  states. 


Threshold  (a:) 

E(TX| 

Transient 

Z0  =  l) 

Asymptotic 

E(TX2 1 

Transient 

Z0  =  l) 

Asymptotic 

0.25 

3.41992E-01 

2.82828E-01 

1.35629E-01 

7.99918E-02 

0.50 

6.45565E-01 

5.65657E-01 

4.74599E-01 

3.19967E-01 

1.00 

1.21845E+00 

1.13131E+00 

1.64830E+00 

1.27987E+00 

5.00 

5.74431E+00 

5.65657E+00 

3.40842E+01 

3.19967E+01 

10.00 

1.14009E+01 

1.13131E+01 

1.32223E+02 

1.27987E+02 

20.00 

2.27140E+01 

2.26263E+01 

5.20481E+02 

5.11948E+02 

50.00 

5.66534E+01 

5.65657E+01 

3.22110E+03 

3.19967E+03 

100.00 

1.13219E+02 

1.13131E+02 

1.28416E+04 

1.27987E+04 

200.00 

2.26350E+02 

2.26263E+02 

5.12807E+04 

5.11948E+04 

500.00 

5.65744E+02 

5.65657E+02 

3.20182E+05 

3.19967E+05 

As  seen  in  Table  1,  the  asymptotic  approximations  track  closely  with  their  transient  counterparts. 
The  discrepancy  between  the  values  is  diminishing  as  x  increases.  In  general,  the  asymptotic  ap¬ 
proximation  underestimates  the  time  to  reach  level  x  and  can,  therefore,  be  viewed  as  a  conservative 
approximation . 

Next  we  compared  the  asymptotic  normal  lifetime  distribution  <h(f;7,  <r)  with  the  numerically 
inverted  lifetime  c.d.f.  using  (33).  Figure  1  displays  three  different  transient  distributions  and 
the  asymptotic  approximation.  For  this  illustration,  we  note  that  the  asymptotic  approximation 
performs  well,  even  when  x  is  small. 
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Figure  1:  Transient  versus  asymptotic  lifetime  c.d.f.s. 

To  further  illustrate  the  quality  of  the  c.d.f.  approximations,  we  compared  the  transient  and 
asymptotic  distributions  using  the  maximum  absolute  deviation  in  probability.  This  error  measure 
is  plotted  against  values  of  x  ranging  from  1.0  to  100.0  in  increments  of  5.0  units  in  Figure  2. 


21 


0.12 


o.i  H 


0.08 


§  0.06 
at 
-a 
a> 


0.02  h 


20 


40 


60 


80 


100 


Figure  2:  Maximum  absolute  deviation  in  probability  as  x  — >  oo. 

Figure  2  indicates  a  possibly  exponential  rate  of  convergence  to  the  limiting  distribution;  however, 
this  cannot  be  verified  analytically.  For  this  example,  the  deviation  is  very  small,  even  with  x  <  10.0. 
Next  we  study  a  more  diverse  environment  with  ten  distinct  states. 

5.2  General  CTMC  Environment 

The  second  example  considers  an  environment  that  transitions  between  ten  states,  i.e.,  the 
environment  {Zt  :  t  >  0}  has  state  space  S  =  {1,  2, . . . ,  10}.  The  10  x  10  infinitesimal  generator 
matrix  is  given  by 


-3.86 

0.06 

0.49 

0.51 

0.07 

0.49 

0.41 

0.93 

0.37 

0.53 

0.39 

-3.37 

0.04 

0.45 

0.38 

0.50 

0.56 

0.26 

0.25 

0.55 

0.25 

0.46 

-3.93 

0.33 

0.37 

0.84 

0.27 

0.20 

0.92 

0.28 

0.35 

0.86 

0.33 

-4.67 

0.48 

0.81 

0.78 

0.05 

0.63 

0.37 

0.74 

0.86 

0.90 

0.89 

-6.18 

0.86 

0.39 

0.61 

0.88 

0.06 

0.65 

0.47 

0.31 

0.76 

0.34 

-4.30 

0.03 

0.55 

0.64 

0.54 

0.94 

0.79 

0.25 

0.88 

0.25 

0.57 

-5.41 

0.10 

0.80 

0.84 

0.83 

0.66 

0.43 

0.46 

0.58 

0.61 

0.56 

-4.71 

0.44 

0.15 

0.47 

0.00 

0.84 

0.80 

0.52 

0.10 

0.20 

0.44 

-3.55 

0.17 

0.63 

0.13 

0.18 

0.13 

0.16 

0.16 

0.09 

0.07 

0.10 

-1.65 
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The  stationary  distribution  of  this  chain  is 


7T  =  [  0.1257  0.1004  0.0894  0.0984  0.0497  0.0992  0.0574  0.0683  0.1179  0.1937  ]• 
The  Markov- modulated  wear  rates  and  shock  arrival  rates  are  summarized  in  Table  2. 


Table  2:  Summary  of  shock  and  wear  rates  for  example  2. 


State  (i) 

Wear  Rate  (r(i)) 

Shock  Arrival  Rate  (A*) 

1 

1.0000 

0.2359 

2 

2.0000 

0.7117 

3 

3.0000 

0.5264 

4 

4.0000 

0.4782 

5 

5.0000 

0.5673 

6 

6.0000 

0.5698 

7 

7.0000 

0.9437 

8 

8.0000 

0.6053 

9 

9.0000 

0.9887 

10 

10.0000 

0.0794 

We  assume  that  the  damage  induced  by  each  shock  follows  an  Erlang  distribution  with  parameters 
2  and  8.0.  That  is,  E(Yi )  =  2/8.0  =  0.25.  The  Laplace-Stieltjes  transform  of  the  shock-magnitude 
c.d.f.  is  therefore 

^(u)  =  (drib)  ’  Re('“)  >  °- 

The  remaining  parameter  values  are  as  follows: 

7  =  1  «  0.16937;  A  =  vrAe  «  0.50731; 

7T  A 

alv  =  13.30623;  erf-  =  0.03125; 

n  =  l/E(Si)  «  0.48523;  a  =  \J +  \cry  ~  2.5411. 

As  before,  we  first  assess  the  quality  of  the  asymptotic  moment  approximations.  The  results  are 
summarized  in  Table  3,  and  they  consistently  show  that  the  asymptotic  approximation  improves 
as  the  degradation  threshold  x  increases.  As  in  the  first  example,  the  asymptotic  approximation 
tends  to  underestimate  the  transient  moments. 
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Table  3:  E(TX  |Zq  =  1),  n  =  1,2,  when  Z  has  ten  states. 


Threshold  (x) 

E(TX| 

Transient 

Z0  =  l) 

Asymptotic 

E(TX2 1 

Transient 

Z0  =  l) 

Asymptotic 

10.00 

1.97116E+00 

1.69365E+00 

4.15994E+00 

2.86846E+00 

25.00 

4.50836E+00 

4.23413E+00 

2.10670E+01 

1.79279E+01 

40.00 

7.03183E+00 

6.77461E+00 

5.07335E+01 

4.58954E+01 

55.00 

9.54684E+00 

9.31509E+00 

9.30839E+01 

8.67709E+01 

70.00 

1.20536E+01 

1.18556E+01 

1.48036E+02 

1.40555E+02 

85.00 

1.45522E+01 

1.43961E+01 

2.15506E+02 

2.07246E+02 

100.00 

1.70425E+01 

1.69365E+01 

2.95412E+02 

2.86846E+02 

115.00 

1.95247E+01 

1.94770E+01 

3.87671E+02 

3.79354E+02 

130.00 

2.19988E+01 

2.20175E+01 

4.92203E+02 

4.84770E+02 

145.00 

2.44647E+01 

2.45580E+01 

6.08924E+02 

6.03094E+02 

160.00 

2.69226E+01 

2.70984E+01 

7.37755E+02 

7.34326E+02 

Next  we  consider  the  lifetime  distribution  approximations.  In  Figure  3  we  plot  three  distinct 
transient  distributions  along  with  the  asymptotic  normal  approximation.  The  discrepancy  between 
the  c.d.f.s  is  significant  when  x  is  small  (namely  when  x  =  1.0),  but  improves  dramatically  as  x 
increases.  When  x  =  100.0,  the  asymptotic  approximation  is  very  accurate  as  expected. 
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Figure  3:  Transient  versus  asymptotic  lifetime  c.d.f.s. 

As  for  the  two-state  example,  we  next  plotted  the  maximum  absolute  deviation  in  probability  as  a 
function  of  the  degradation  threshold,  x.  Consistent  with  the  two-state  case,  Figure  4  shows  that 
the  worst-case  deviation  is  monotonically  decreasing  (to  zero)  as  x  — >  oo. 
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Figure  4:  Transient  versus  asymptotic  lifetime  c.d.f.s. 

The  two  numerical  examples  illustrate  the  usefulness  of  the  asymptotic  approximations,  namely 
that  they  allow  one  to  circumvent  the  task  of  numerically  inverting  Laplace  transforms.  The  next 
section  summarizes  our  main  contributions  and  provides  some  concluding  remarks. 

6  Conclusions 

In  this  paper  we  have  presented  both  transient  and  asymptotic  reliability  indices  for  a  single¬ 
unit  system  whose  degradation  process  may  be  modeled  as  a  Markov-modulated  shock  and  wear 
process.  This  analytical  framework  provides  a  great  deal  of  modeling  flexibility,  particularly  in  the 
context  of  degradation-based  reliability  modeling,  by  allowing  for  a  wide  range  of  wear  dynamics 
and  the  inclusion  of  damage-inducing  shocks  that  are  dictated  by  the  device’s  random  operating 
environment.  These  characteristics  make  it  appealing  for  modeling  real  systems  when  degradation 
can  be  directly  observed,  or  when  a  direct  mapping  between  the  environmental  conditions  and  the 
degradation  can  be  accurately  modeled  (e.g.,  by  employing  specific  physics-of- failure  models). 

While  the  models  are  mathematically  sound  and  easy  to  implement,  in  their  current  form  they 
lack  the  flexibility  to  account  for  environment  state  sojourn  times,  or  shock  inter-arrival  times,  that 
are  not  exponentially  distributed.  In  this  sense,  the  models  are  somewhat  restrictive.  Moreover, 
the  models  assume  that  future  wear  and  damage  are  independent  of  the  history  of  the  degradation 
process  -  an  assumption  that  may  be  very  difficult  to  justify  in  practice.  Nevertheless,  the  main 
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results  provide  a  framework  for  simple,  asymptotically  accurate  approximations  that  can  serve  as 
a  starting  point  for  applications  requiring  real-time  updating  of  remaining  lifetime  distributions 
(as  in  a  condition-based  maintenance  environment).  Specifically,  real-time  observations  of  the 
evolution  of  the  environment,  or  the  degradation  level,  can  be  used  to  update  the  parameters 
of  the  modulating  process,  the  state-dependent  wear  rates,  and  the  shock  arrival  rates,  using 
the  techniques  described  in  Kharoufeh  and  Cox  [14].  Tracking  the  occurrence  of  environment 
transitions,  and  the  number  of  shocks  occurring  over  a  (sufficiently  long)  time  interval  will  provide 
reasonable  statistical  estimates,  so  long  as  the  number  of  environment  states  is  not  too  large. 
The  existence  of  a  minimal  representation  of  the  modulating  Markov  chain  that  leads  to  accurate 
distribution  estimates  was  posited  in  [14].  Such  a  representation  indicates  that  the  results  of  this 
work  can  be  implemented  in  a  realistic  setting. 

In  the  future  it  will  be  instructive  to  investigate  similar  approximation  schemes  that  allow 
for  a  generalized  environment  process  that  does  not  evolve  as  a  continuous-time  Markov  chain. 
If  the  models  are  to  be  of  any  practical  value,  it  will  be  necessary  to  accurately  estimate  the 
functions  r  and  A  describing  the  evolution  of  wear  and  damage  as  a  function  of  the  environment. 
For  this  purpose,  real  degradation  data  is  required,  as  is  the  guidance  and  experience  of  subject 
matter  experts,  to  ensure  the  wear  and  damage  models  are  appropriate  for  the  application.  An 
examination  of  multi-unit  systems  is  also  warranted.  However,  for  a  multi-unit  system  in  a  random 
environment,  it  may  not  be  possible  to  treat  each  unit  independently  because,  although  the  unit 
lifetimes  may  be  statistically  independent,  they  share  a  common  dependence  on  the  operating 
environment;  therefore,  under  certain  conditions,  the  unit  lifetimes  are  associated  to  one  another 
as  shown  by  Qinlar  et  al.  [4] .  Nevertheless,  our  approach  allows  for  a  multi- unit  system  to  be  viewed 
as  a  single-unit  system  if  the  aggregate  influence  of  the  environment  can  be  discerned.  Finally,  an 
important  issue  that  we  did  not  consider  in  this  paper  is  the  rate  of  convergence  of  the  lifetime 
moments  and  distributions  to  their  respective  limits.  Our  empirical  results  indicate  that  the  rate 
of  convergence  may  be  exponential,  though  this  conjecture  cannot  be  confirmed  analytically  using 
the  transform  results.  It  will  be  instructive  to  derive  error  bounds  for  our  approximations  in  future 
work. 
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